Libraries

library(survival)
library(KMsurv)
library(ggplot2)
library(ggpubr)
library(survminer)
library(plotly)
library(muhaz)
library(ggthemes)

Loading in Data

data(burn)

burn1 <- burn
burn1 <-
  data.frame(burn1, Treatment = factor(burn1$Z1, labels = c("Routine", "Cleansing")))
burn1 <-
  data.frame(burn1, Gender = factor(burn1$Z2, labels = c("Male", "Female")))
burn1 <-
  data.frame(burn1, Race = factor(burn1$Z3, labels = c("Nonwhite", "White")))
burn1 <- data.frame(burn1, PercentBurned = burn1$Z4)
burn1 <-
  data.frame(burn1, SiteHead = factor(burn1$Z5, labels = c("NotBurned", "Burned")))
burn1 <-
  data.frame(burn1, SiteButtock = factor(burn1$Z6, labels = c("NotBurned", "Burned")))
burn1 <-
  data.frame(burn1, SiteTrunk = factor(burn1$Z7, labels = c("NotBurned", "Burned")))
burn1 <-
  data.frame(burn1, SiteUpperLeg = factor(burn1$Z8, labels = c("NotBurned", "Burned")))
burn1 <-
  data.frame(burn1, SiteLowerLeg = factor(burn1$Z9, labels = c("NotBurned", "Burned")))
burn1 <-
  data.frame(burn1, SiteRespTract = factor(burn1$Z10, labels = c("NotBurned", "Burned")))
burn1 <-
  data.frame(burn1, BurnType = factor(burn1$Z11, labels = c("Chemical", "Scald", "Electric", "Flame")))

burn1.surv <- with(burn1, Surv(T3, D3))
# plot(survfit(burn1.surv ~ Treatment, data = burn1),
#      col = 1:2,
#      lwd = 2)
# title("Time to Infection for Routine Care and Total Body Cleansing")
# legend(
#   "topright",
#   c("Routine Care", "Total Body Cleansing"),
#   col = 1:2,
#   lwd = 2
# )

#Lets observe the data 
burn1
print(survdiff(burn1.surv ~ Treatment, data = burn1))
## Call:
## survdiff(formula = burn1.surv ~ Treatment, data = burn1)
## 
##                      N Observed Expected (O-E)^2/E (O-E)^2/V
## Treatment=Routine   70       28     21.4      2.07      3.79
## Treatment=Cleansing 84       20     26.6      1.66      3.79
## 
##  Chisq= 3.8  on 1 degrees of freedom, p= 0.05

Question 1

# Plot the Kaplan Meier Curves 
KMcurves <- survfit(burn1.surv ~ Treatment, data = burn1)

KMplot <-
  ggsurvplot(KMcurves, ggtheme = theme_fivethirtyeight()) + labs(title = "KM Curves for Time to Infection for Routine Care \n and Total Body Cleansing")
ggplotly(KMplot[[1]])

Question 2

cumHazPlot <-
  ggsurvplot(
    KMcurves,
    fun = "cumhaz",
    conf.int = FALSE,
    palette = c("#5d8aa8", "#e32636"),
    ggtheme =theme_fivethirtyeight()
  ) + ggtitle("Cumulative Hazard for Treatment Type")

#ggplotly(cumHazPlot[[1]])

#cumHazPlot
ggplotly(cumHazPlot$plot) # + geom_abline())
## Warning: plotly.js does not (yet) support horizontal legend items 
## You can track progress here: 
## https://github.com/plotly/plotly.js/issues/53
#Cloglog plot
cLogLogPlot <-
  ggsurvplot(
    KMcurves,
    fun = "cloglog",
    palette = c("#5d8aa8", "#e32636"),
    ggtheme = theme_fivethirtyeight()
  ) + ggtitle("Complimentary Log-Log for Treatment") 
ggplotly(cLogLogPlot[[1]])
## Warning: plotly.js does not (yet) support horizontal legend items 
## You can track progress here: 
## https://github.com/plotly/plotly.js/issues/53

Question 3

#Building with time independent objects
survTimeIndep <- Surv(burn1$T3, burn1$D3)
#Build a cox model with just treatment
burn.cox.indep.T <- coxph(survTimeIndep ~ Treatment, data = burn1)
summary(burn.cox.indep.T)
## Call:
## coxph(formula = survTimeIndep ~ Treatment, data = burn1)
## 
##   n= 154, number of events= 48 
## 
##                       coef exp(coef) se(coef)      z Pr(>|z|)  
## TreatmentCleansing -0.5614    0.5704   0.2934 -1.914   0.0557 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    exp(coef) exp(-coef) lower .95 upper .95
## TreatmentCleansing    0.5704      1.753     0.321     1.014
## 
## Concordance= 0.566  (se = 0.039 )
## Likelihood ratio test= 3.73  on 1 df,   p=0.05
## Wald test            = 3.66  on 1 df,   p=0.06
## Score (logrank) test = 3.76  on 1 df,   p=0.05
#Building with treatement and burnpercentage
# burn.cox.indep.TPB <- coxph(survTimeIndep ~ Treatment + PercentBurned, data = burn1)
# summary(burn.cox.indep.TPB)

#Building with treatement, burnpercentage, and burn Type 
# burn.cox.indep.TPBbt <- coxph(survTimeIndep ~ Treatment + PercentBurned + BurnType, data = burn1)
# summary(burn.cox.indep.TPBbt)


#Building with treatement, burnpercentage, and burn Type + Respartory Tract
# burn.cox.indep.TPBbtResp <- coxph(survTimeIndep ~ Treatment + PercentBurned + BurnType + SiteRespTract, data = burn1)
# summary(burn.cox.indep.TPBbtResp)

#Building with treatement, burnpercentage, and burn Type + resp Tract and gender
# burn.cox.indep.TPBbtRespG <- coxph(survTimeIndep ~ Treatment + PercentBurned + BurnType + SiteRespTract + Gender, data = burn1)
# summary(burn.cox.indep.TPBbtRespG)

# Build a model with all of the above + race
summary(coxph(survTimeIndep ~ Treatment + PercentBurned + BurnType + SiteRespTract + Gender + Race, data = burn1))
## Call:
## coxph(formula = survTimeIndep ~ Treatment + PercentBurned + BurnType + 
##     SiteRespTract + Gender + Race, data = burn1)
## 
##   n= 154, number of events= 48 
## 
##                          coef exp(coef)  se(coef)      z Pr(>|z|)  
## TreatmentCleansing  -0.622888  0.536393  0.303140 -2.055   0.0399 *
## PercentBurned        0.002904  1.002909  0.007712  0.377   0.7065  
## BurnTypeScald        1.582955  4.869324  1.099898  1.439   0.1501  
## BurnTypeElectric     2.090905  8.092234  1.089444  1.919   0.0550 .
## BurnTypeFlame        0.915903  2.499032  1.024396  0.894   0.3713  
## SiteRespTractBurned  0.216496  1.241718  0.359962  0.601   0.5475  
## GenderFemale        -0.559528  0.571479  0.402071 -1.392   0.1640  
## RaceWhite            2.250113  9.488803  1.031151  2.182   0.0291 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                     exp(coef) exp(-coef) lower .95 upper .95
## TreatmentCleansing     0.5364     1.8643    0.2961    0.9717
## PercentBurned          1.0029     0.9971    0.9879    1.0182
## BurnTypeScald          4.8693     0.2054    0.5639   42.0440
## BurnTypeElectric       8.0922     0.1236    0.9566   68.4549
## BurnTypeFlame          2.4990     0.4002    0.3356   18.6097
## SiteRespTractBurned    1.2417     0.8053    0.6132    2.5143
## GenderFemale           0.5715     1.7498    0.2599    1.2567
## RaceWhite              9.4888     0.1054    1.2575   71.6026
## 
## Concordance= 0.738  (se = 0.039 )
## Likelihood ratio test= 24.8  on 8 df,   p=0.002
## Wald test            = 19.9  on 8 df,   p=0.01
## Score (logrank) test = 23.44  on 8 df,   p=0.003
# We can see from p-values that Treatment, Burn Type, and Race are the only variables with significance, so we will build a cox model using these
burn.cox <- coxph(survTimeIndep ~ Treatment + BurnType + Race, data = burn1)
summary(burn.cox)
## Call:
## coxph(formula = survTimeIndep ~ Treatment + BurnType + Race, 
##     data = burn1)
## 
##   n= 154, number of events= 48 
## 
##                       coef exp(coef) se(coef)      z Pr(>|z|)  
## TreatmentCleansing -0.6009    0.5483   0.2979 -2.018   0.0436 *
## BurnTypeScald       1.5785    4.8477   1.0865  1.453   0.1463  
## BurnTypeElectric    2.1719    8.7745   1.0856  2.001   0.0454 *
## BurnTypeFlame       1.0035    2.7278   1.0163  0.987   0.3234  
## RaceWhite           2.2842    9.8182   1.0260  2.226   0.0260 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    exp(coef) exp(-coef) lower .95 upper .95
## TreatmentCleansing    0.5483     1.8238    0.3058     0.983
## BurnTypeScald         4.8477     0.2063    0.5764    40.775
## BurnTypeElectric      8.7745     0.1140    1.0450    73.676
## BurnTypeFlame         2.7278     0.3666    0.3722    19.992
## RaceWhite             9.8182     0.1019    1.3144    73.340
## 
## Concordance= 0.705  (se = 0.037 )
## Likelihood ratio test= 21.85  on 5 df,   p=6e-04
## Wald test            = 17.4  on 5 df,   p=0.004
## Score (logrank) test = 20.61  on 5 df,   p=0.001

LRT and Wald Test both provide evidence that there is a significant difference with patients with differrent treatments, burn types, and race.

Question 4

# We check Residuals

#Schoenfield
burn.zph <- cox.zph(burn.cox)


ggcoxzph(
  burn.zph[1],
  se = FALSE,
  font.main = 12,
  ggtheme = theme_fivethirtyeight(),
  point.col = "#5d8aa8",
)

ggcoxzph(
  burn.zph[2],
  se = FALSE,
  font.main = 12,
  ggtheme = theme_fivethirtyeight(),
  point.col = "#5d8aa8",
)

ggcoxzph(
  burn.zph[3],
  se = FALSE,
  font.main = 12,
  ggtheme = theme_fivethirtyeight(),
  point.col = "#5d8aa8"
)

# Martingale Residuals 
burn.mart<- residuals(burn.cox, type = "martingale")

burn.lowess <- as.data.frame(lowess(burn1$T3, burn.mart))

# Plot Martingale 
ggplotly(
  ggplot() + aes(burn1$T3, burn.mart) + geom_point(col = "#5d8aa8") + labs(x = "Time", y = "Martingale Residuals", title = "Martingale Residuals vs. Time till Infection with Straphylocous Aureaus") + geom_line(data = burn.lowess, aes(x = x, y = y), col = "#388E3C") + theme_fivethirtyeight()
)
# We can see that we need a transformation, so we use log time 
burn.lowess.log <-  as.data.frame(lowess(log(burn1$T3), burn.mart))

# Plot Martingale 
ggplotly(
  ggplot() + aes(log(burn1$T3), burn.mart) + geom_point(col = "#5d8aa8") + labs(x = "Time", y = "Martingale Residuals", title = "Martingale Residuals vs. Time till Infection with Straphylocous Aureaus") + geom_line(data = burn.lowess.log, aes(x = x, y = y), col = "#388E3C") + theme_fivethirtyeight()
)
# Let us make time categorical 
burn1["CatTime"] <-
  as.factor(cut(burn1$T3,
      c(0, 49, 100),
      labels = c("<50", "50+")))

#Buld a new cox Model with new time variable 
burn.cox.catTime <- coxph(survTimeIndep ~ Treatment + BurnType + Race + CatTime, data = burn1)
## Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
## Loglik converged before variable 6 ; coefficient may be infinite.
summary(burn.cox.catTime)
## Call:
## coxph(formula = survTimeIndep ~ Treatment + BurnType + Race + 
##     CatTime, data = burn1)
## 
##   n= 154, number of events= 48 
## 
##                          coef  exp(coef)   se(coef)      z Pr(>|z|)  
## TreatmentCleansing -6.226e-01  5.366e-01  2.966e-01 -2.099   0.0358 *
## BurnTypeScald       1.918e+00  6.809e+00  1.087e+00  1.765   0.0775 .
## BurnTypeElectric    2.200e+00  9.023e+00  1.086e+00  2.025   0.0428 *
## BurnTypeFlame       1.099e+00  3.000e+00  1.017e+00  1.080   0.2800  
## RaceWhite           2.539e+00  1.266e+01  1.031e+00  2.462   0.0138 *
## CatTime50+         -1.906e+01  5.271e-09  3.880e+03 -0.005   0.9961  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    exp(coef) exp(-coef) lower .95 upper .95
## TreatmentCleansing 5.366e-01  1.864e+00    0.3000    0.9596
## BurnTypeScald      6.809e+00  1.469e-01    0.8095   57.2769
## BurnTypeElectric   9.023e+00  1.108e-01    1.0734   75.8500
## BurnTypeFlame      3.000e+00  3.333e-01    0.4088   22.0154
## RaceWhite          1.266e+01  7.897e-02    1.6775   95.5835
## CatTime50+         5.271e-09  1.897e+08    0.0000       Inf
## 
## Concordance= 0.743  (se = 0.035 )
## Likelihood ratio test= 37.73  on 6 df,   p=1e-06
## Wald test            = 18.65  on 6 df,   p=0.005
## Score (logrank) test = 28.58  on 6 df,   p=7e-05
### DOES NOT CONVERGE

Question 5

# Now we add time dependent covariates
burn.cox.TD <- coxph(survTimeIndep ~ Treatment + BurnType + Race + T1 + D1+T2+D2 , data = burn1)
summary(burn.cox.TD)
## Call:
## coxph(formula = survTimeIndep ~ Treatment + BurnType + Race + 
##     T1 + D1 + T2 + D2, data = burn1)
## 
##   n= 154, number of events= 48 
## 
##                         coef exp(coef)  se(coef)      z Pr(>|z|)  
## TreatmentCleansing -0.351516  0.703621  0.314807 -1.117   0.2642  
## BurnTypeScald       1.440387  4.222331  1.097976  1.312   0.1896  
## BurnTypeElectric    1.783797  5.952413  1.111258  1.605   0.1084  
## BurnTypeFlame       0.831245  2.296176  1.021213  0.814   0.4157  
## RaceWhite           2.304509 10.019254  1.029397  2.239   0.0252 *
## T1                  0.001372  1.001373  0.021945  0.063   0.9502  
## D1                 -0.351020  0.703970  0.437191 -0.803   0.4220  
## T2                  0.001054  1.001054  0.015185  0.069   0.9447  
## D2                 -0.828456  0.436723  0.457322 -1.812   0.0701 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    exp(coef) exp(-coef) lower .95 upper .95
## TreatmentCleansing    0.7036    1.42122    0.3796     1.304
## BurnTypeScald         4.2223    0.23684    0.4909    36.320
## BurnTypeElectric      5.9524    0.16800    0.6742    52.553
## BurnTypeFlame         2.2962    0.43551    0.3103    16.993
## RaceWhite            10.0193    0.09981    1.3323    75.346
## T1                    1.0014    0.99863    0.9592     1.045
## D1                    0.7040    1.42052    0.2988     1.658
## T2                    1.0011    0.99895    0.9717     1.031
## D2                    0.4367    2.28978    0.1782     1.070
## 
## Concordance= 0.754  (se = 0.036 )
## Likelihood ratio test= 30.51  on 9 df,   p=4e-04
## Wald test            = 25.62  on 9 df,   p=0.002
## Score (logrank) test = 29.51  on 9 df,   p=5e-04
# We remove the ones that are not significant to the model
# Now we add time dependent covariates
burn.cox.TD <- coxph(survTimeIndep ~ Treatment + BurnType + Race + D2 , data = burn1)
summary(burn.cox.TD)
## Call:
## coxph(formula = survTimeIndep ~ Treatment + BurnType + Race + 
##     D2, data = burn1)
## 
##   n= 154, number of events= 48 
## 
##                       coef exp(coef) se(coef)      z Pr(>|z|)   
## TreatmentCleansing -0.4218    0.6559   0.3034 -1.390  0.16455   
## BurnTypeScald       1.5112    4.5324   1.0884  1.388  0.16500   
## BurnTypeElectric    2.0049    7.4252   1.0867  1.845  0.06505 . 
## BurnTypeFlame       0.8808    2.4129   1.0173  0.866  0.38658   
## RaceWhite           2.3147   10.1214   1.0282  2.251  0.02438 * 
## D2                 -0.9012    0.4061   0.3386 -2.662  0.00777 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    exp(coef) exp(-coef) lower .95 upper .95
## TreatmentCleansing    0.6559     1.5246    0.3619    1.1888
## BurnTypeScald         4.5324     0.2206    0.5368   38.2649
## BurnTypeElectric      7.4252     0.1347    0.8825   62.4771
## BurnTypeFlame         2.4129     0.4144    0.3285   17.7206
## RaceWhite            10.1214     0.0988    1.3490   75.9383
## D2                    0.4061     2.4625    0.2091    0.7885
## 
## Concordance= 0.751  (se = 0.036 )
## Likelihood ratio test= 29.57  on 6 df,   p=5e-05
## Wald test            = 24.37  on 6 df,   p=4e-04
## Score (logrank) test = 28.06  on 6 df,   p=9e-05

Question 6

drop1(burn.cox.TD, test = "Chisq")
# We check Residuals

#Schoenfield
burn.zph.TD <- cox.zph(burn.cox.TD)


ggcoxzph(
  burn.zph.TD[1],
  se = FALSE,
  font.main = 12,
  ggtheme = theme_fivethirtyeight(),
  point.col = "#5d8aa8",
)

ggcoxzph(
  burn.zph.TD[2],
  se = FALSE,
  font.main = 12,
  ggtheme = theme_fivethirtyeight(),
  point.col = "#5d8aa8",
)

ggcoxzph(
  burn.zph.TD[3],
  se = FALSE,
  font.main = 12,
  ggtheme = theme_fivethirtyeight(),
  point.col = "#5d8aa8"
)

ggcoxzph(
  burn.zph.TD[4],
  se = FALSE,
  font.main = 12,
  ggtheme = theme_fivethirtyeight(),
  point.col = "#5d8aa8"
)